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Case  Study  of  Model-Based 
Inversion  of  the  Angle  Beam 
Ultrasonic  Response  From 
Composite  Impact  Damage 

The  U.S.  Air  Force  seeks  to  improve  lifecycle  management  of  composite  structures.  Non¬ 
destructive  characterization  of  damage  is  a  key  input  to  this  framework.  One  approach  to 
characterization  is  model-based  inversion  of  ultrasound  inspection  data;  however,  the 
computational  expense  of  simulating  the  response  from  damage  represents  a  major  hur¬ 
dle  for  practicality.  A  surrogate  forward  model  with  greater  computational  efficiency  and 
sufficient  accuracy  is,  therefore,  critical  to  enable  damage  characterization  via  model- 
based  inversion.  In  this  work,  a  surrogate  model  based  on  Gaussian  process  regression 
(GPR)  is  developed  on  the  chirplet  decomposition  of  the  simulated  quasi-shear  scatter 
from  delamination-like  features  that  form  a  shadowed  region  within  a  representative 
composite  layup.  The  surrogate  model  is  called  in  the  solution  of  the  inverse  problem  for 
the  position  of  the  hidden  delamination,  which  is  achieved  with  <0.5%  error  in  <20  min 
on  a  workstation  computer  for  two  unique  test  cases.  These  results  demonstrate  that  solv¬ 
ing  the  inverse  problem  from  the  ultrasonic  response  is  tractable  for  composite  impact 
damage  with  hidden  delaminations.  [DOT:  TO. 1115/1.4040233] 


1  Introduction 

The  increasing  use  of  polymer-matrix  composites  within  Air 
Force  structures  has  fueled  a  growing  demand  for  improved  life- 
cycle  management  of  these  materials  [1],  The  Air  Force  Aircraft 
Structural  Integrity  Program  guidelines  for  metallics  provide  the 
template  for  this  improvement,  where  validated  damage  progres¬ 
sion  models  enable  tolerance  to  assumed  defects  for  an  acceptable 
measure  of  risk  [2],  Composites  damage  progression  models  are 
sensitive  to  the  geometric  complexity  of  typical  damage  modes. 
Consequently,  nondestructive  evaluation  (NDE)  techniques  are 
needed  to  volumetrically  represent  damage  for  improved  lifecycle 
management. 

Currently  viable  composite  NDE  methods  include  X-ray- 
computed  tomography  (XCT)  and  ultrasound.  XCT  scans  volu¬ 
metrically  reconstruct  internal  damage,  making  it  an  obvious  path 
for  feature  characterization.  Nevertheless,  XCT  cannot  be  used  in 
situ  as  a  single-sided  inspection  and  generally  requires  additional 
safety  precautions.  By  contrast,  normal  incidence  longitudinal 
wave  ultrasound  can  detect  the  presence  of  delaminations  and 
matrix  cracking  in  a  single-sided  inspection  (Fig.  1(a)).  However, 
it  is  incapable  of  insonifying  the  hidden  damaged  region  under  the 
top  delamination  layers  (Fig.  1(b)).  A  novel  angle  beam  shear 
wave  method  addresses  this  limitation  by  using  direct  and  back- 
wall  reflected  quasi-shear  waves  to  inspect  for  hidden  ply  delami¬ 
nations  and  matrix  cracks  [3,4].  This  approach  builds  on  prior 
work  on  polar  backscatter  ultrasonic  techniques  [5-8].  In  this 
way,  additional  information  from  hidden  damage  features  can  be 
captured  and  processed  to  characterize  the  damage  producing  the 
observed  response.  It  should  be  noted  that  phased  array  trans¬ 
ducers  (either  in  a  normal  incidence  or  angle  beam  shear  configu¬ 
ration)  may  provide  additional  information;  however,  for  the 
purposes  of  this  work,  the  simplest  sensing  approach  (single 
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element,  pulse-echo  inspection)  was  selected  based  on  the  com¬ 
plexity  of  the  material  and  internal  damage. 

One  challenge  with  this  angle  beam  ultrasonic  inspection  tech¬ 
nique  is  resolving  the  multiple  signals  scattered  from  delamination 
edges  at  different  depths.  Manual  interpretation  of  the  paths  for 
each  signal  is  complicated,  as  the  scattered  signals  interact  with 
the  composite  panel  surfaces  and  neighboring  delaminations. 
Model-based  inversion  is  proposed  to  address  this  challenge.  This 
approach  uses  nonlinear  least-squares  based  optimization  to  esti¬ 
mate  the  parameters  of  a  physics-based  forward  model  by  mini¬ 
mizing  an  objective  function  defined  as  the  squared  residual 
between  the  model  output  and  the  observed  response.  Prior  work 
has  considered  flaw  reconstructions  with  ultrasonic  NDE  in  aniso¬ 
tropic  materials  [9-13],  Other  efforts  have  examined  the  use  of 
multipath  signals  for  characterizing  the  location  of  scatterers 
[14-17],  For  example,  Hutt  and  Simonetti  used  a  strongly  scatter¬ 
ing  planar  interface  as  a  mirror  to  look  behind  the  object  to 
achieve  complete  reconstruction  [17].  Some  studies  have  investi¬ 
gated  the  identification  of  multiply  diffracted  ultrasonic  waves 
from  point  scatterers,  although  many  of  them  have  focused  on 
medical  ultrasound  applications  [18,19].  Most  approaches  out¬ 
lined  within  the  literature  have  leveraged  approximate  models  like 
ray  theory  for  improving  reconstruction.  The  proposed  solution 
builds  on  prior  work  using  numerical  solutions  of  parameterized 
models  to  invert  the  damage  state  [20-22].  The  specific  problem 
of  inverting  the  location  of  multiple  impact  damage  discontinu¬ 
ities  in  a  thin  composite  panel  has  not  been  addressed  to  date.  In 
this  work,  a  simplified  case  study  was  investigated  to  demonstrate 
the  model-based  inversion  approach. 


2  Materials  and  Methods 

The  geometry  of  impact  damage  within  a  composite  is  depend¬ 
ent  on  the  size,  shape,  and  velocity  of  the  impactor;  the  interaction 
of  the  deforming  impactor  with  the  shape  and  orientation  of  the 
structure;  and  the  current  mechanical  state  of  the  material  as  a 
function  of  as-manufactured  properties,  internal  defects,  thermo¬ 
mechanical  degradation,  and  environment.  These  variables  result 
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in  asymmetric,  highly  irregular  damage  contours.  Bench-top 
impact  tests  on  laboratory  coupons  can  control  some  of  these  vari¬ 
ables,  albeit  with  little  reduction  in  the  geometric  complexity  of 
the  damage. 

A  simple  representative  problem,  estimation  of  the  location  of  a 
hidden  delamination  beneath  a  top  delamination,  was  considered 
to  reduce  geometric  complexity  to  a  level  suitable  for  a  feasibility 
study.  Two  cases  are  presented  in  Fig.  2:  in  case  1,  the  lower  (hid¬ 
den)  delamination  was  fixed  in  z  with  the  x-position  iterated 
through  nine  values;  in  case  2,  the  lower  delamination  was  fixed 
in  x  with  the  z-position  iterated  through  eight  values. 

The  parameters  for  the  study  follow  existing  specimens  and  lab¬ 
oratory  testing  conditions  [23].  A  3.2  mm  thick,  24  ply  IM7/977-3 
composite  panel  of  [—45  deg/90  deg/45  deg/0  deg]3s  layup  with 
simplified  delamination-like  features  was  considered.  The  simu¬ 
lated  ultrasonic  transducer  was  6.35  mm  in  diameter,  with  a  center 
frequency  of  5  MHz  and  a  focal  length  of  19.05  mm.  The  angle  of 
incidence  was  set  to  24  deg,  near  the  first  critical  angle  of  the 
homogenized  composite  (~27  deg)  such  that  the  incident  beam 
was  refracted  into  quasi-shear  modes  at  the  front-wall.  The  simu¬ 
lated  delaminations  were  each  12.7  mm  long.  Their  locations  did 
not  directly  coincide  with  ply  interfaces  but  did  not  exceed  a  mini¬ 
mum  of  a  single  ply  separation  in  the  z-direction.  These  cases  pro¬ 
vide  an  elegant  demonstration  of  the  proposed  method  of  solving 
the  inverse  problem  while  leaving  room  for  follow-on  studies 


Fig.  1  Impact  delamination  in  a  composite  coupon.  Amplitude 
and  time-of-flight  data  describes  the  complexity  of  composite 
impact  damage  and  hidden  delamination  regions  invisible  to  a 
normal  incidence  longitudinal  wave  single-sided  inspection:  (a) 
normal  incidence  C-Scan  of  impact  damage  with  visible  petal¬ 
shaped  deiaminations  based  on  amplitude  data  and  (b)  repre¬ 
sentation  of  a  hidden  region  (white)  in  a  delamination  field 
(black)  based  on  time-of-flight  data. 
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addressing  the  challenges  of  a  much  larger  set  of  pertinent 
parameters. 

3  Calculation 

3.1  Forward  Model  and  Need  for  Surrogate  Representa¬ 
tion.  The  efficiency  of  model-based  inversion  is  proportional  to 
the  well-posedness  of  the  problem  and  inversely  proportional  to 
the  cost  of  each  forward  model  iterate.  Two  computational 
tools — CIVA  FIDEL  2D  and  PZFlex — are  tailored  for  ultrasound 
simulation  of  composite  materials  [24-26].  CIVA  FIDEL  2D  cou¬ 
ples  semi-analytical  beam  calculations  with  a  finite  difference 
time  domain  (FDTD)  solver  acting  on  a  user-specified  region  to 
model  two-dimensional  (2D)  wave  propagation  for  simple  geome¬ 
tries  and  boundary  conditions  [27,28],  Comparatively,  PZFlex 
uses  finite  element  analysis  with  an  optimized  explicit  solver  to 
simulate  2D  or  three-dimensional  wave  propagation  on  any  geom¬ 
etry  with  complex  boundary  conditions.  Both  tools  can  simulate  a 
spatially  dense  B-scan  response  from  thin  delamination-like  fea¬ 
tures  in  a  representative  composite  material  within  several  hours; 
however,  typical  inverse  problems  often  require  hundreds  of 
model  calls.  For  a  2D  ultrasonic  inspection  scenario,  solution 
times  would  range  from  days  to  weeks.  When  the  physics-based 
model  is  computationally  expensive,  a  surrogate  will  often  be 
developed  that  sacrifices  some  accuracy  for  a  decrease  in 


X 


1  ft  nam. 

(a) 


(b) 


Fig.  2  Schematics  of  surrogate  model  demonstration  cases  1 
(a)  and  2  (b).  The  transducer  is  6.35  mm  in  diameter,  with  a  cen¬ 
ter  frequency  of  5  MHz  and  a  19.05  mm  focal  length.  Delamina¬ 
tion  features  are  12.7  mm  long.  The  diagramed  scan  path  is  not 
to  scale:  (a)  hidden  delamination  varying  in  x-direction,  fixed  z 
at  z=2mm  and  (b)  hidden  delamination  varying  in  z-direction, 
fixed  xat  x=  1.2  mm. 
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(a)  Transducer  position  (mm)  (^)  Transducer  position  (mm)  (c)  Transducer  position  (mm) 

Fig.  3  Chirplet  representation  of  an  example  B-scan.  (a)  The  original  B-scan  from  CIVA.  (b)  The  chirplet  reconstruction  of 
the  B-scan.  (c)  The  residual  between  the  original  B-scan  and  the  reconstruction.  The  color  bar  represents  the  amplitude  of 
the  response  in  (a)  and  (b)  and  the  amplitude  of  the  residual  in  (c)  [34]. 


computational  expense  [29].  The  need  for  an  efficient  surrogate 
forward  model  is  evident  for  the  ultrasound  inverse  problem;  thus, 
CIVA  was  chosen  for  surrogate  model  development  based  on  its 
greater  computational  efficiency,  the  dimensionality  of  the  sce¬ 
nario,  and  previous  validation  work  [4]. 


Previous  work  demonstrated  the  feasibility  of  an  ultrasound 
surrogate  model  for  the  case  of  a  single  delamination  placed  in 
four  unique  z-positions  [30],  In  this  work,  a  similar  model  is 
developed  on  the  more  challenging  case  of  a  hidden  delamination. 
Model  development  began  with  selection  of  the  representative 


( g )  Transducer  Position  (mm)  (/))  Transducer  Position  (mm)  (/)  Transducer  Position  (mm) 

Fig.  4  Hidden  delamination  simulations  for  nine  hidden  delamination  x-positions.  The  x-position  range  was  selected  to  ensure 
that  the  delamination  was  hidden  but  not  so  far  from  the  tip  of  the  topmost  delamination  that  indications  from  the  hidden 
delamination  were  no  longer  visible.  The  color  bar  describes  the  amplitude  of  the  response:  (a)  x=  0.8  mm,  (b)  x=  1.0  mm,  (c) 
x=  1.2  mm,  (d)  x=  1.4  mm,  (e)  x=  1.6  mm,  (f)  x  =  1.8  mm,  (g)  x=  2.0  mm,  (/?)  x=  2.2  mm,  and  (#)  x=  2.4  mm. 
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Fig.  5  Gaussian  chirplet  parameters  and  GPR  estimates  for  case  1.  Open  circles  are  parameters  from  CIVA  and 
closed  circles  are  parameter  estimates,  (a)  and  (c)  vary  smoothly,  while  (e)  and  (f)  exhibit  strong  variation,  (a) 
Reflection  1  amplitude,  (b)  reflection  2  bandwidth,  (c)  reflection  1  time  of  arrival,  (d)  reflection  2  center  frequency, 
(e)  reflection  8  phase,  and  (f)  reflection  3  chirp  rate. 
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problem,  which  was  modeled  in  CIV  A  FIDEL  2D  for  a  sparse  set 
of  design  points  containing  the  inverse  solution.  Chirplet  decom¬ 
position  was  applied  to  each  simulated  B-scan;  then,  the  surrogate 
model  was  constructed  by  interpolating  on  the  chirplet  parameter 
space  via  Gaussian  process  regression  (GPR).  The  resulting  fits 
formed  the  surrogate  model  used  to  solve  the  inverse  problem  for 
the  ultrasonic  response  from  a  hidden  delamination  with  an 
unknown  x-  or  z-position. 

3.2  B-Scan  Simulation.  Several  assumptions  were  made  for 
the  simulation.  First,  it  was  assumed  that  the  composite  was  flat 
and  parallel,  which  is  valid  for  a  simulation-based  study  but  likely 
invalid  for  real  composite  panels.  Second,  it  was  assumed  that 
delaminations  could  be  modeled  as  discontinuities  with  zero 
thickness  parallel  to  the  plies.  While  acceptable  for  a  simulation- 
based  study,  real  delamination  morphologies  should  be  incorpo¬ 
rated  to  validate  this  assumption.  Finally,  it  was  assumed  that 
water  is  present  on  the  backside  of  the  composite.  This  assump¬ 
tion  would  hold  true  for  a  real  panel  scanned  via  immersion  ultra¬ 
sound;  however,  an  air  back-wall  is  more  likely  for  real  inspection 
scenarios. 

The  model  consisted  of  a  water  path  to  simulate  an  immersion- 
based  inspection,  the  plies  forming  the  composite  panel,  and 
delamination  features.  Because  of  the  FIDEL  package,  the  plies 
were  modeled  individually  (the  composite  was  not  homogenized). 
For  these  boundary  conditions,  a  2D  FDTD  region  was  set  to  sur¬ 
round  the  near  edges  of  the  delaminations  and  the  composite  panel 
boundaries.  The  FDTD  solution  region  moves  with  the  scanned 
probe,  so  the  lateral  extent  was  set  to  a  large  size  (20  mm)  in  order 
to  include  both  direct  and  full  skip  interactions  with  the  far-wall. 
Perfectly  matched  layers  were  also  applied  to  the  lateral  bounda¬ 
ries  to  help  control  fictitious  signals  from  the  model  domain 
edges.  A  mesh  resolution  was  fixed  to  l/20th  of  the  center  fre¬ 
quency  wavelength.  Additional  details  of  the  formulation  can  be 
found  in  published  literature  [27,28].  B-scan  images  were 


generated  by  scanning  the  transducer  over  a  6  mm  path  across  the 
top  of  the  composite  in  0. 1  mm  increments  with  a  constant  water 
path  of  17  mm.  Simulations  were  performed  on  a  HP  computer 
with  two  12-core  hyperthreaded  processers  (48  compute  cores)  at 
2.5  GHz  clock  speed  and  98  GB  of  RAM.  Each  B-scan  required 
~360min  of  simulation  time.  An  example  simulated  B-scan  for  a 
single  delamination  located  at  1.0  mm  in  depth  from  the  top  sur¬ 
face  is  shown  in  Fig.  3(fl).  In  this  case,  the  largest  signal  is  a  direct 
diffraction  from  the  delamination  edge.  Signals  later  in  time  corre¬ 
spond  to  the  multiple  paths  for  diffracted  responses  interacting 
with  the  top  or  bottom  surfaces  of  the  composite.  Half-skips 
involve  a  single  reflection  off  of  the  back-wall  surface,  while  full- 
skips  involve  two  reflections  off  of  the  back-wall.  These  signals 
shift  in  time  and  space  according  to  the  position  of  the  delamina¬ 
tion  edge. 

3.3  Chirplet  Decomposition.  It  is  nontrivial  to  produce  B- 
scans  via  the  physics-based  forward  model  given  the  nature  of  the 
shifting  transient  signals  with  varying  delamination  position  and 
the  size  of  the  simulated  B-scans  results  (61  scan  steps  x  1350 
time  steps  ~  80k  data  points).  To  reduce  the  dimensionality  of  for¬ 
ward  model  results,  individual  A-scans  were  represented  as  a  lin¬ 
ear  combination  of  Gaussian  chirplets  [31].  The  Gaussian  chirplet 
is  defined  as 

/(f)  =  jSexp^— oq  ( t  —  t)~  jcos^27t/c(t  —  t)  +  (f>  +  1x2(1  —  t)“^ 

(1) 

where  /  is  the  amplitude  of  the  Gaussian  envelope,  oq  is  the  band¬ 
width  of  the  Gaussian  envelope,  t  is  the  time  of  arrival  of  the 
response,  fc  is  the  center  frequency,  cj>  is  the  phase  angle,  and  oc2  is 
the  chitp  rate.  The  chitp  rate,  a2>  captures  the  variation  of  fre¬ 
quency  with  respect  to  time  and  represents  a  substantial  improve¬ 
ment  over  a  wavelet-based  approach.  The  chirplet  has  been  shown 
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Fig.  6  Comparison  of  CIVA  and  surrogate  model  B-scans:  (a)-(c)  x=  1.2  mm  and  ( d)-(f)  x=2.0mm.  (a)  and  ( d)  depict  the 
original  CIVA  B-scan,  ( b )  and  (e)  depict  the  chirplet  reconstruction  of  the  B-scan,  and  (c)  and  (/)  depict  the  residual  between 
the  CIVA  B-scan  and  the  chirplet  reconstruction.  The  residual  is  a  mixture  of  untracked  reflections  and  incomplete  chirplet 
fitting  and  is  everywhere  low.  (a)  CIVA  B-scan,  x=  1.2  mm,  (£>)  chirplet  reconstruction,  x=  1.2  mm,  (c)  residual,  x=  1.2  mm,  (d) 
CIVA  B-scan,  x=  2.0  mm,  (e)  chirplet  reconstruction,  x=  2.0  mm,  and  (/)  residual,  x=  2.0  mm. 
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to  model  individual  ultrasonic  reflections  as  a  sparse,  energy  pre¬ 
serving  representation  of  the  data  without  any  a  priori  information 
on  the  expected  response  [31].  The  chirplet  parameters  can  be 
estimated  from  the  analytic  signal  and  further  refined  using  a 
Gauss-Newton  algorithm  [32].  In  the  chirplet  decomposition 
algorithm  suggested  in  Refs.  [31-33],  chirplets  are  fit  successively 
to  individual  A-scans  until  a  reconstruction  error  tolerance  is 
reached.  The  method  used  here  differs  slightly;  rather  than  faith¬ 
fully  reconstructing  each  A-scan,  chirplets  only  reproduced  reflec¬ 
tions  in  the  B-scan  that  represented  responses  from  a 
delamination.  This  was  achieved  by  screening  out  reflections  from 
the  front-wall  and  back-wall  of  the  composite — the  remaining 
reflections  were  visually  segmented  and  fitted.  The  surrogate 
model  predicts  how  these  individual  reflections  change  with  the 
position  of  the  hidden  delamination.  Chirplet  decomposition 
reduces  the  number  of  parameters  needed  to  represent  the  B-scan 
to  rij,  where  N  is  the  total  number  of  distinct  reflections, 

and  rij  is  the  number  of  transducer  positions  containing  a  signal 
from  the  y'th  reflection.  As  an  example.  Fig.  3  depicts  a  B-scan 
containing  ~80k  samples,  while  the  chirplet  reconstruction  is 
defined  by  ~360  chirplet  parameters  [34]. 

3.4  Gaussian  Process  Regression.  Techniques  such  as  splin- 
ing,  support  vector  machines,  and  GPR  were  considered  for  use  as 
a  surrogate  for  the  physics-based  forward  model.  Gaussian  pro¬ 
cess  regression,  also  known  as  kriging,  was  chosen  as  the  meta¬ 
model  for  its  efficiency,  robustness,  and  flexibility,  and  has  a 
history  of  use  as  a  metamodel  [35-37].  A  type  of  supervised 


learning,  GPR  generates  input-output  mappings  from  the  training 
data.  The  estimation  of  the  function  value  /  at  a  test  point  u  is 
described  by 

f(u)  =  ji(u)  +  cl  (C  +  (fy)  \y-n{s))  (2) 

The  assumed  mean  function,  /t,  describes  the  general  trend  of  the 
data.  Both  the  covariance  vector  cu  between  the  test  point(s)  u  and 
the  sample  points  ,v  and  the  covariance  matrix  C  between  the  sam¬ 
ple  points  .v  depend  on  the  assumed  covariance  function.  In  the 
previous  formulation,  the  known  responses  y  recorded  at  the  sam¬ 
ple  points  s  are  assumed  to  be  the  sum  of  the  unknown  latent  func¬ 
tion  values  f(s)  and  independent,  identically  distributed  Gaussian 
noise  with  variance  an2.  For  this  application,  the  inputs  .?  are  given 
by  the  transducer  location  and  x-  or  z-position  of  the  delamination, 
while  the  outputs  y  are  the  value  of  the  chirplet  parameter  /,  oq,  t, 
fc,  (j),  or  a2.  For  the  case  of  simulation  data,  the  variance  a  2  is 
assumed  to  be  zero.  The  outputs /are  the  estimates  of  the  chirplet 
parameters  for  the  transducer  locations  at  intermediate  x-  or  z- 
positions  of  the  hidden  delamination. 

Development  of  the  GPR  model  begins  with  selection  of  the 
mean  and  covariance  function.  Both  are  assumed  to  be  stationary. 
Example  functions  tested  for  fitting  are  presented  for  the  mean 
(Eq.  (3))  and  covariance  (Eq.  (4)) 

{  a0  +  a\u  +  L  +  a„u"  }  ^ 


Fig.  7  CIVA  B-scan  (a),  inverse  solution  ( b ),  surrogate  model  evaluated  at  the  inverse  solution  (c),  and 
the  residual  in  decibels  (d).  The  inverse  solution  has  an  error  of  0.22%,  and  the  surrogate  model  at  the 
inverse  solution  compares  very  well  with  the  CIVA  B-scan.  (a)  CIVA  B-scan  for  x  =  1.3  mm,  ( b )  inverse 
solution  for  50  initial  guesses  over  20  iterations,  (c)  surrogate  model  evaluated  at  x  =  1 .297  mm,  and  (d) 
residual  in  decibels. 
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k(Si,Sj) 


exp(-||j;  -  Sy||2//) 
exp^—  ||s,-  —  s;||2/2/2j 


(4) 


In  Eq.  (3),  u  are  the  test  point(s),  c  is  a  constant,  a,  are  polynomial 
coefficients,  and  n  is  the  order  of  the  term.  In  Eq.  (4),  s  are  the 
sample  points  and  ||  ||2  denotes  the  Euclidean  distance.  The 
parameter  /  describes  how  rapidly  the  correlation  between  samples 
decreases  and  is  a  hyperparameter  of  the  GPR  model  optimized 
via  leave-one-out  cross  validation  on  the  input  data;  thus,  the 
value  of  /  varies  between  parameters  and  reflections  [38].  Selec¬ 
tion  of  the  mean  and  covariance  functions  is  often  guided  by 
physics-based  expectations  of  the  response;  however,  mathemati¬ 
cally  rigorous  model  selection  techniques  such  as  the  Akaike 
information  criterion  (AIC,  Eq.  (5))  can  be  applied 

AIC  =  — (2/A)£[log  life]  +  (2  d/N)  (5) 

The  AIC  rewards  model  accuracy  (the  first  term,  where  N  is  the 
number  of  samples)  while  penalizing  complexity  (the  second 
term,  where  d  is  the  number  of  model  parameters),  with  better 


models  producing  more  highly  negative  values  [39].  For  this 
work,  the  response  from  each  parameter  was  fit  with  every  combi¬ 
nation  of  the  selected  mean  and  covariance  functions,  with  the 
combination  producing  the  lowest  AIC  selected  as  the  optimal  set 
for  that  parameter.  The  surrogate  model  results  from  the  best  per¬ 
forming  mean  and  covariance  function  combination  with  an  opti¬ 
mized  hyperparameter  I.  An  input  vector  u  is  developed  for  each 
intermediate  x-  or  z-position  of  the  hidden  delamination  and 
applied  to  the  GPR  model.  The  output  of  the  model  is  the  estimate 
of  the  selected  parameter  for  a  given  reflection. 


4  Results  and  Discussion 

4.1  Case  1:  x-Varying,  z-Fixed.  For  this  demonstration,  the 
z-position  of  the  lower  delamination  was  fixed  at  2  mm  from  the 
upper  delamination,  while  x  ranged  from  0.8  mm  to  2.4  mm  in 
nine  steps  from  the  leftmost  edge  of  the  upper  delamination.  This 
range  was  selected  due  to  the  evident  variations  in  reflection  posi¬ 
tion  and  intensity  at  each  delamination  position.  Changes  between 
B-scans  diminish  drastically  beyond  x  =  2.4  mm,  resulting  in 
regions  where  unique  solutions  to  the  inverse  problem 


Fig.  8  Hidden  deiamination  simulations  for  eight  second-delamination  z-positions.  The  z-position  range  was  selected  to  ensure 
that  the  second  delamination  was  hidden  but  not  so  close  to  the  topmost  delamination  that  the  reflections  fully  merge  and 
become  indistinguishable.  The  color  bar  describes  the  amplitude  of  the  response:  (a)  z=  1.1  mm,  (b)  z=  1.3  mm,  (c)  z=  1.5  mm, 
(d)  z=  1.7  mm,  (e)  z=  1.9  mm,  (/)  z=  2.1  mm,  (g)  z  =  2.3  mm,  and  ( h )  z  =  2.5  mm. 
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hypothetically  do  not  exist.  The  temporospatial  translations  of 
nine  reflections  were  tracked  across  each  B-scan.  These  reflections 
were  selected  in  order  of  decreasing  amplitude,  with  the  bounds 
of  each  reflection  governed  by  the  envelope  of  the  chirplet  in  the 
time  dimension  and  the  threshold  for  fitting  a  chirplet  to  an 
observed  signal  in  the  spatial  direction.  The  B-scans  used  for 
building  the  surrogate  model  are  presented  in  Fig.  4.  Tracked 
reflections  generally  maintain  their  original  shape  as  the  hidden 
delamination  varies  from  x=0.8mm  to  x=  1.4mm,  although  the 
strength  of  the  response  from  later  reflections  begins  to  diminish. 
From  x=  1.4  to  x=  2.0  more  reflections  appear,  while  those  later 
in  time  continue  to  decrease  in  amplitude.  Finally,  later  reflections 
diminish  to  the  point  where  their  amplitude  is  below  the  threshold 
for  chirplet  fitting  between  x  =  2.0  and  x  =  2.4. 

Sample  data  were  drawn  from  Gaussian  chirplet  parameters  for 
a  given  delamination  x-position.  Three  mean  basis  functions  (con¬ 
stant,  linear,  and  quadratic)  (Eq.  (3))  and  two  covariance  basis 
functions  (Eq.  (4))  were  competed  for  each  parameter  fit  using  the 
AIC.  In  general,  the  Gaussian  covariance  function  provided  the 
most  accurate  representation  of  the  observed  data.  The  optimal 
mean  function  varied  from  parameter  to  parameter,  although  less 
sensitivity  to  the  mean  function  was  observed.  To  demonstrate  the 
surrogate  model  for  case  1 ,  selected  Gaussian  chirplet  parameters 
and  the  GPR  model  fits  from  the  nine  reflections  for  x  =  1.3  mm 
are  presented  in  Fig.  5.  Several  of  the  parameters  vary  smoothly 
across  a  given  reflection  (e.g.,  amplitude  and  time  of  arrival,  Figs. 
5(a)  and  5(c),  respectively),  although  few  parameters  are  uni¬ 
formly  smooth  across  all  reflections  (e.g.,  reflection  eight  phase 
and  reflection  three  chirp  rate,  Figs.  5(c)  and  5(f),  respectively). 
As  the  quality  of  the  chirplet  representation  is  most  sensitive  to 
the  amplitude  and  time  of  arrival  of  the  signal,  the  smoothness  of 
the  response  from  these  parameters  is  highly  encouraging  for 
accurate  representation  of  an  intermediate  B-scan  location. 

Fitted  parameters  were  used  to  reconstruct  the  B-scan  for  test 
positions  of  x=  1.2  mm  and  r  =  2.0mm  as  presented  in  Figs. 


6(a)-6(c)  and  6(d)-6(f),  respectively.  The  quality  of  the  surrogate 
was  judged  by  the  residual  between  the  simulated  B-scan  for  that 
delamination  location  and  the  estimated  B-scan.  Construction  of 
the  B-scan  from  the  surrogate  model  took  52  ms.  As  the  nine 
tracked  reflections  were  selected  based  in  order  of  decreasing 
amplitude,  the  residual  for  both  cases  is  dominated  by  low- 
amplitude  reflections;  however,  some  response  from  tracked 
reflections  (most  notably  the  residual  arising  from  the  first,  largest 
reflection)  was  not  fully  modeled  by  the  chirplet  fits. 

The  true  performance  of  the  surrogate  model  is  measured  by  its 
ability  to  rapidly  and  accurately  invert  a  B-scan  with  a  delamina¬ 
tion  in  an  unknown  position.  This  was  tested  by  first  simulating  a 
CIV  A  B-scan  for  an  intermediate  point  in  the  surrogate  model 
(lower  delamination  positioned  atxm  1.3  mm).  Next,  an  optimiza¬ 
tion  algorithm  was  applied  to  minimize  the  residual  sum-of- 
squares  between  the  supplied  B-scan  and  that  generated  by  the 
surrogate  model.  As  reflections  may  appear  or  disappear  as  the  x- 
position  of  the  hidden  delamination  is  changed,  an  algorithm  that 
handles  discontinuities  in  the  search  space  was  required.  The  dif¬ 
ferential  evolution  (DE)  algorithm  was,  thus,  selected  [40]. 

The  CIVA  B-scan,  initial  guesses  and  final  solution  of  the  DE 
algorithm,  surrogate  model  evaluated  at  the  final  solution,  and  the 
residual  in  decibels  are  presented  in  Figs.  l(a)-l(d),  respectively. 

Twenty  iterations  of  the  DE  algorithm  were  evaluated  using  50 
initial  points  covering  the  possible  solution  space,  resulting  in  a 
final  solution  of  x  =  1.297  mm.  This  result  has  an  error  of  0.22% 
with  a  runtime  of  ~  18  min  on  a  standard  workstation.  Visually, 
the  B-scan  reconstructed  for  x=  1.297  mm  appears  nearly  indis¬ 
tinguishable  from  the  CIVA  B-scan  at  x=  1.3  mm,  further  demon¬ 
strating  the  quality  of  the  inversion. 

4.2  Case  2:  z-Varying,  x-Fixed.  For  case  2,  the  x-position  of 
the  lower  delamination  was  fixed  at  1.2  mm  from  the  leftmost 
edge  of  the  upper  delamination,  while  z  ranged  from  1 . 1  mm  to 


((/j  Transducer  Position  (mm)  (q'j  Transducer  Position  (mm)  Transducer  Position  (mm) 

Fig.  9  Comparison  of  CIVA  and  surrogate  model  B-scans,  (a)-(c)  z  =  1.3  mm  and  (d)-(f)  z=  2.1  mm.  (a)  and  ( d)  depict  the  origi¬ 
nal  CIVA  B-scan,  ( b )  and  (e)  depict  the  chirplet  reconstruction  of  the  B-scan,  and  (c)  and  (f)  depict  the  residual  between  the  CIVA 
B-scan  and  the  chirplet  reconstruction.  The  residual  is  a  mixture  of  untracked  reflections  and  incomplete  chirplet  fitting  and  is 
everywhere  low:  (a)  CIVA  B-scan,  z=  1.3  mm,  (b)  chirplet  reconstruction,  z=  1.3  mm,  (c)  residual,  z=  1.3  mm,  (cf)  CIVA  B-scan, 
z=  2.1  mm,  (e)  chirplet  reconstruction,  z=  2.1  mm,  and  (f)  residual,  z=  2.1  mm. 
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Fig.  10  CIVA  B-scan  (a),  inverse  solution  ( b ),  surrogate  model  evaluated  at  the  inverse  solution  ( c ),  and  the  resid¬ 
ual  in  decibels  (d).  The  inverse  solution  has  an  error  of  0.41%,  and  the  surrogate  model  at  the  inverse  solution  com¬ 
pares  very  well  with  the  CIVA  B-scan.  (a)  CIVA  B-scan  for  z=  2.2  mm,  (b)  inverse  solution  for  50  initial  guesses  over 
20  iterations,  (c)  surrogate  model  evaluated  at  z=  2.191  mm,  and  ( d)  residual  in  decibels. 


2.5  mm  from  the  upper  delamination  in  eight  steps.  As  before,  the 
range  of  hidden  delamination  z-positions  was  chosen  to  ensure 
distinct  variations  in  the  intensity  and  position  of  the  tracked 
reflections.  The  temporospatial  translations  of  nine  reflections 
were  tracked  across  each  B-scan.  Reflections  were  selected  in 
order  of  decreasing  amplitude,  with  the  bounds  of  each  reflection 
governed  by  the  envelope  of  the  chirplet  in  the  time  dimension 
and  the  threshold  for  fitting  a  chirplet  to  an  observed  signal  in  the 
spatial  direction.  The  B-scans  developed  for  the  surrogate  model 
are  presented  in  Fig.  8.  As  the  hidden  delamination  moves  from 
zm  1.1  mm  to  z=  1.5  mm,  the  reflections  around  27  fis  increase  in 
size  and  complexity,  although  the  general  shape  of  the  response 
remains  mostly  constant.  Complex  reflections  beyond  25  /ts 
appear  and  shift  with  the  --position  of  the  hidden  delamination 
from  z  =  1 .5  mm  to  z  =  2. 1  mm. 

Fitted  parameters  were  used  to  reconstruct  the  B-scan  for  test 
positions  of  z=  1.3  mm  and  z  =  2.1mm  as  presented  in  Figs. 
9(a)-9(c)  and  9(d)-9(f),  respectively.  The  quality  of  the  surrogate 
was  judged  by  the  residual  between  the  simulated  B-scan  for  that 
delamination  location  and  the  estimated  B-scan. 

Construction  of  the  B-scan  from  the  surrogate  model  took 
61  ms.  As  in  case  1,  the  residual  is  dominated  by  low  amplitude 
reflections  or  incomplete  fitting  of  the  primary  reflection. 

To  test  this  case,  a  CIVA  B-scan  was  generated  for  the  lower 
delamination  positioned  at  z  =  2.2mm.  The  CIVA  B-scan,  initial 
guesses  and  final  solution  of  the  DE  algorithm,  surrogate  model 
evaluated  at  the  inverse  solution,  and  residual  in  decibels  are  pre¬ 
sented  in  Figs.  10(a)-10(<f),  respectively.  The  DE  algorithm  was 


again  used  with  20  iterations  on  50  initial  points  covering  the  pos¬ 
sible  solution  space.  An  inverse  solution  of  z  =  2.191  was  found. 
This  result  has  error  of  0.41%  with  a  runtime  of  ~14min  on  a 
standard  workstation.  As  observed  in  case  1,  the  B-scan  recon¬ 
structed  for  z  =  2.191  mm  appears  nearly  indistinguishable  from 
the  CIVA  B-scan  at  z  =  2.2  mm,  further  demonstrating  the  quality 
of  the  inversion. 

5  Conclusions 

The  ultrasonic  response  from  a  composite  material  with  embed¬ 
ded  delamination-like  features  was  simulated  for  two  test  cases. 
Gaussian  chirplets  were  fit  to  selected  reflections  within  the  simu¬ 
lated  B-scans.  Chirplet  parameters  and  the  position  of  the  shifting 
coordinate  ( x  or  z)  were  used  to  develop  GPR  surrogate  models. 
These  models  were  used  to  estimate  the  unknown  position  of  the 
hidden  delamination  for  each  test  case.  The  unknown  positions 
were  rapidly  estimated  to  within  <0.5%  of  the  actual  value,  dem¬ 
onstrating  that  an  efficient  surrogate  model  can  be  developed  on 
ultrasound  simulations  and  that  inversion  of  the  ultrasonic 
response  from  delamination-like  features  is  tractable  for  compos¬ 
ite  impact  damage  with  hidden  delaminations. 

Practical  application  begins  with  an  initial  model  of  the  mate¬ 
rial  instantiated  with  the  general  shape  of  the  delamination  field 
based  on  a  priori  knowledge  of  relevant  parameters  (impact 
energy,  layup,  etc.).  B-scans  for  permutations  of  this  field  would 
be  generated,  decomposed  by  the  chirplet  transform,  and  parame¬ 
terized  with  Kriging  to  produce  a  surrogate  model.  The  residual 
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between  the  parameterized  model  and  the  B-scan  of  the  actual 
damage  would  be  minimized,  resulting  in  damage  features  match¬ 
ing  those  within  the  real  composite  structure. 

Several  challenges  drive  future  work.  First,  chirplet  decomposi¬ 
tion  relies  on  accurately  mapping  distinct  reflections  within  the  B- 
scan.  While  some  degree  of  separation  between  responses  is 
imposed  by  layup  and  ply  thickness,  the  proximity  of  damage  fea¬ 
tures  often  results  in  merged  reflections.  The  current  mapping 
approach  requires  user  interpretation  of  the  total  response,  which 
only  compounds  the  challenge  of  selecting  reflections  to  fit  within 
the  surrogate  model.  Algorithmic  identification  of  individual 
reflections  is,  thus,  crucial  to  enabling  the  inverse  solution  for 
complex  damage  morphologies.  Second,  the  inverse  problem  was 
demonstrated  for  design  points  where  reflections  from  the  hidden 
delamination  remain  within  the  B-scan.  These  indications  vanish 
as  the  hidden  feature  becomes  further  occluded  by  the  upper 
delamination.  Successive  positions  of  the  hidden  delamination 
will  eventually  result  in  the  same  observed  response,  resulting  in 
nonunique  solutions  to  the  inverse  problem.  Better  methods  of 
insonifying  the  hidden  region  is  thereby  critical  to  extending  the 
range  of  tractable  design  cases.  Third,  experimental  noise  and  var¬ 
iability  from  factors  such  as  probe  misalignment,  water  path,  and 
transducer  representation  (e.g.,  dimensions/properties  of  the  mod¬ 
eled  piezoelectric  crystal,  backing  material,  or  lens)  are  antici¬ 
pated  to  impact  the  quality  of  the  inverse  solution.  Parameter 
sensitivity  studies  will  be  employed  to  understand  the  impact  of 
these  sources  of  uncertainty.  Finally,  the  techniques  developed 
within  this  work  could  be  extended  to  other  relevant  composite 
damage  features,  including  fiber  breakage  and  matrix  cracking. 
Each  feature  type  dominates  failure  for  specific  composite  config¬ 
urations,  rendering  them  equally  important  in  predicting  the 
remaining  loading  cycles  survivable  by  the  damaged  component. 
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